use ${newdata}carownership_dataset_bergen, clear

/* Only keep required obs */
keep if $trmgroup != 0 
keep if $trmgroup != .
keep if 	(year >= $firstyearpre  & year <= $lastyearpre) ///
		  | (year >= $firstyearpost & year <= $lastyearpost)

/* setup variables */
capt drop bergen 		  
gen commuters = ($trmgroup == 1 | $trmgroup == 3) 
gen bergen = 	($trmgroup == 1 | $trmgroup == 2)
gen post = (year >= $firstyearpost & year <= $lastyearpost)
gen bergen_commuters = bergen * commuters
gen commuters_post = commuters * post
gen bergen_commuters_post = bergen_commuters * post

/* All time-varying variables based on geography (residence or work location)
   are set to their 2014 values (variation post 2014 is endogenous to treatment) */
foreach var in dist time_work PublicVSCarTime_fam_mean PublicDiffCarTime_fam_mean grk {
	replace `var' = . if year != 2014
	bysort familienr (`var'): replace `var' = `var'[_n-1] if missing(`var')
}

sort familienr year
egen long famid_num = group(familienr)

/* === First dimension - income groups. Total household income ============== */
gen heterovar_inc = 0
replace heterovar_inc = wies1 if wies1 != .
replace heterovar_inc = heterovar_inc + wies2 if wies2 != .	

* Creating two income groups. top/bottom 20 %
xtile inc_qt2 = heterovar_inc if year == 2014, nquantiles(5)
replace inc_qt2 = . if inc_qt2 != 1 & inc_qt2 != 5
replace inc_qt2 = 2 if inc_qt2 == 5	

/* === High/low income, quintiles of public transport measure =============== */		  

* Groups: 1-5 based on public transport vs car time to work
xtile heterovar_ptdiff_qt = PublicDiffCarTime_fam_mean if year == 2014, nquantiles(5)
* Groups: 1-10 by interacting with one of two income groups
gen ptdiff_inc_qt2 = heterovar_ptdiff_qt*2 - 2 + inc_qt2
	
bys fam_fe: egen aux = mean(ptdiff_inc_qt2)
replace ptdiff_inc_qt2 = aux
drop aux*

/* === Running regression ====================================================*/
capt drop bergen_commuters_post_*	
forvalues i = 1/10 {
	gen bergen_commuters_post_`i' = bergen_commuters_post * (ptdiff_inc_qt2 == `i')
}			  

eststo clear
reghdfe bev ///
	bergen_commuters_post_* /// DiDiD: post commuter in bergen
	i.year##i.ptdiff_inc_qt2 /// 
	$xvar ///
	i.commuters##i.year##i.ptdiff_inc_qt2 /// 
	i.bergen_commuters##i.ptdiff_inc_qt2 ///
	i.bergen##i.ptdiff_inc_qt2 ///
	, absorb( ///
		i.grk#i.year ///
		famid_num ///
	) vce(cluster grk)	

/* === Storing results in variables ==========================================*/	
matrix B = e(b)
matrix V = e(V)

capt drop xaxis* beta1 beta2 se1 se2
gen xaxis = .
gen beta1 = .
gen se1 = .
gen beta2 = .
gen se2 = .

forvalues i = 1/5 {
	local i1 = `i'* 2 - 1
	local i2 = `i'* 2
	qui replace xaxis = `i' in `i'
	qui replace beta1  = B[1,`i1'] in `i'
	qui replace beta2  = B[1,`i2'] in `i'
	qui replace se1    = sqrt(V[`i1',`i1']) in `i'
	qui replace se2    = sqrt(V[`i2',`i2']) in `i'
}

/* 95 % CI */
capt drop plusse* minusse*
qui gen plusse1 = beta1 + 1.96 * se1
qui gen minusse1 = beta1 - 1.96 * se1
qui gen plusse2 = beta2 + 1.96 * se2
qui gen minusse2 = beta2 - 1.96 * se2
gen xaxis2 = xaxis + 0.1
gen xaxis1 = xaxis - 0.1


/* === FIGURE ============================================================*/

twoway ///	
	(rcap plusse1 minusse1 xaxis1 ///
		, sort color(blue)) ///			
	(scatter beta1 xaxis1 ///
		, sort lcolor(black) mcolor(blue) lpattern(solid) msymbol(O)) ///	
	(rcap plusse2 minusse2 xaxis2 ///
		, sort color(green)) ///			
	(scatter beta2 xaxis2 ///
		, sort lcolor(black) mcolor(green) lpattern(solid) msymbol(D)) ///		
	, graphregion(color(white))	///
	xtitle("Quintiles of public transit quality") ytitle("Treatment effect") ///
	name(figE1c, replace) ///
	xlabel(1 "< 29 min" 2 "29-37 min" 3 "37-48 min" 4 "48-71 min" 5 "> 71 min", angle(55)) ///
	legend(on order(2 4) ///
		label(2 "Low income") ///
		label(4 "High income") ///
		pos(6) col(2) ///
	) ///
	yline(0) ylabel(-0.04(0.04)0.16, angle(0)) scale(1.3)
	
graph export "${figures}figE1c.png", as(png) replace
graph export "${figures}figE1c.pdf", as(pdf) replace		
		

	
/* === High low income, levels of education ================================= */		  

gen heterovar_educ = .
replace heterovar_educ = max_educ if year == 2014
replace heterovar_educ = . if heterovar_educ == 0
gen educ_inc_qt2 = heterovar_educ*2 - 2 + inc_qt2
bys famid_num: egen aux = mean(educ_inc_qt2)
replace educ_inc_qt2 = aux
drop aux*

/* === Running regression ====================================================*/
capt drop bergen_commuters_post_*	
forvalues i = 1/10 {
	gen bergen_commuters_post_`i' = bergen_commuters_post * (educ_inc_qt2 == `i')
}			  

eststo clear
reghdfe bev ///
	bergen_commuters_post_* /// DiDiD: post commuter in bergen
	i.year##i.educ_inc_qt2 /// 
	$xvar ///
	i.commuters##i.year##i.educ_inc_qt2 /// 
	i.bergen_commuters##i.educ_inc_qt2 ///
	i.bergen##i.educ_inc_qt2 ///
	, absorb( ///
		i.grk#i.year ///
		famid_num ///
	) vce(cluster grk)	

/* === Storing results in variables ==========================================*/	
matrix B = e(b)
matrix V = e(V)

capt drop xaxis* beta1 beta2 se1 se2 
gen xaxis = .
gen beta1 = .
gen se1 = .
gen beta2 = .
gen se2 = .

forvalues i = 1/4 {
	local i1 = `i'* 2 - 1
	local i2 = `i'* 2
	qui replace xaxis = `i' in `i'
	qui replace beta1  = B[1,`i1'] in `i'
	qui replace beta2  = B[1,`i2'] in `i'
	qui replace se1    = sqrt(V[`i1',`i1']) in `i'
	qui replace se2    = sqrt(V[`i2',`i2']) in `i'
}

/* 95 % CI */
capt drop plusse* minusse*
qui gen plusse1 = beta1 + 1.96 * se1
qui gen minusse1 = beta1 - 1.96 * se1
qui gen plusse2 = beta2 + 1.96 * se2
qui gen minusse2 = beta2 - 1.96 * se2
gen xaxis2 = xaxis + 0.1
gen xaxis1 = xaxis - 0.1

/* === FIGURE ============================================================*/

twoway ///	
	(rcap plusse1 minusse1 xaxis1 ///
		, sort color(blue)) ///			
	(scatter beta1 xaxis1 ///
		, sort lcolor(black) mcolor(blue) lpattern(solid) msymbol(O)) ///	
	(rcap plusse2 minusse2 xaxis2 ///
		, sort color(green)) ///			
	(scatter beta2 xaxis2 ///
		, sort lcolor(black) mcolor(green) lpattern(solid) msymbol(D)) ///		
	, graphregion(color(white))	///
	/*xtitle("Quintiles of income")*/ ytitle("Treatment effect") ///
	name(figE1a, replace) ///
	xlabel(1 "< high school" 2 "High school" 3 "Uni., short" 4 "Uni., long", angle(55)) ///
	legend(on order(2 4) ///
		label(2 "Low income") ///
		label(4 "High income") ///		
		pos(6) col(2) ///
	) ///
	/*legend(off)*/ yline(0) ylabel(-0.04(0.04)0.16, angle(0)) scale(1.3)
	
graph export "${figures}figE1a.png", as(png) replace
graph export "${figures}figE1a.pdf", as(pdf) replace	

/* === high/low income, levels of work distance ============================= */		  

xtile dist_qt = dist if year == 2014, nquantiles(5)
gen dist_inc2 = dist_qt*2 - 2 + inc_qt2
bys famid_num: egen aux = mean(dist_inc2)
replace dist_inc2 = aux
drop aux*

* === Run regression ===========================================================

capt drop bergen_commuters_post_*	
forvalues i = 1/10 {
	gen bergen_commuters_post_`i' = bergen_commuters_post * (dist_inc2 == `i')
}			  

eststo clear
reghdfe bev ///
	bergen_commuters_post_* /// DiDiD: post commuter in bergen
	i.year##i.dist_inc2 /// 
	$xvar ///
	i.commuters##i.year##i.dist_inc2 /// 
	i.bergen_commuters##i.dist_inc2 ///
	i.bergen##i.dist_inc2 ///
	, absorb( ///
		i.grk#i.year ///
		famid_num ///
	) vce(cluster grk)	
	
/* === Storing results in variables ==========================================*/

matrix B = e(b)
matrix V = e(V)

/* Storing them as variables */
capt drop xaxis* beta1 beta2 se1 se2 
gen xaxis = .
gen beta1 = .
gen se1 = .
gen beta2 = .
gen se2 = .

forvalues i = 1/5 {
	local i1 = `i'* 2 - 1
	local i2 = `i'* 2
	qui replace xaxis = `i' in `i'
	qui replace beta1  = B[1,`i1'] in `i'
	qui replace beta2  = B[1,`i2'] in `i'
	qui replace se1    = sqrt(V[`i1',`i1']) in `i'
	qui replace se2    = sqrt(V[`i2',`i2']) in `i'
}

/* 95 % CI */
capt drop plusse* minusse*
qui gen plusse1 = beta1 + 1.96 * se1
qui gen minusse1 = beta1 - 1.96 * se1
qui gen plusse2 = beta2 + 1.96 * se2
qui gen minusse2 = beta2 - 1.96 * se2
gen xaxis2 = xaxis + 0.1
gen xaxis1 = xaxis - 0.1

/* === FIGURE ================================================================*/


twoway ///	
	(rcap plusse1 minusse1 xaxis1 ///
		, sort color(blue)) ///			
	(scatter beta1 xaxis1 ///
		, sort lcolor(black) mcolor(blue) lpattern(solid) msymbol(O)) ///	
	(rcap plusse2 minusse2 xaxis2 ///
		, sort color(green)) ///			
	(scatter beta2 xaxis2 ///
		, sort lcolor(black) mcolor(green) lpattern(solid) msymbol(D)) ///		
	, graphregion(color(white))	///
	xtitle("Quintiles of work distance") ytitle("Treatment effect") ///
	name(figE1d, replace) ///
	xlabel(1 "5-7 km" 2 "7-10 km" 3 "10-13 km" 4 "13-19 km" 5 "19-50 km", angle(55)) ///
	legend(on order(2 4) ///
		label(2 "Low income") ///
		label(4 "High income") ///
		pos(6) col(2) ///
	) ///
	/*legend(off)*/ yline(0) ylabel(-0.04(0.04)0.16, angle(0)) scale(1.3)
	
graph export "${figures}figE1d.png", as(png) replace
graph export "${figures}figE1d.pdf", as(pdf) replace

/* === High/low income, quintiles of age ==================================== */		  


xtile age_qt = age if year == 2014, nquantiles(5)
gen age_inc2 = age_qt*2 - 2 + inc_qt2
bys fam_fe: egen aux = mean(age_inc2)
replace age_inc2 = aux
drop aux*

/* === Run regression ======================================================= */

capt drop bergen_commuters_post_*	
forvalues i = 1/10 {
	gen bergen_commuters_post_`i' = bergen_commuters_post * (age_inc2 == `i')
}			  

eststo clear
reghdfe bev ///
	bergen_commuters_post_* /// DiDiD: post commuter in bergen
	i.year##i.age_inc2 /// 
	$xvar ///
	i.commuters##i.year##i.age_inc2 /// 
	i.bergen_commuters##i.age_inc2 ///
	i.bergen##i.age_inc2 ///
	, absorb( ///
		i.grk#i.year ///
		famid_num ///
	) vce(cluster grk)	

/* === Storing results as variables ========================================= */	
matrix B = e(b)
matrix V = e(V)

capt drop xaxis* beta1 beta2 se1 se2 
gen xaxis = .
gen beta1 = .
gen se1 = .
gen beta2 = .
gen se2 = .

forvalues i = 1/5 {
	local i1 = `i'* 2 - 1
	local i2 = `i'* 2
	qui replace xaxis = `i' in `i'
	qui replace beta1  = B[1,`i1'] in `i'
	qui replace beta2  = B[1,`i2'] in `i'
	qui replace se1    = sqrt(V[`i1',`i1']) in `i'
	qui replace se2    = sqrt(V[`i2',`i2']) in `i'
}

/* 95 % CI */
capt drop plusse* minusse*
qui gen plusse1 = beta1 + 1.96 * se1
qui gen minusse1 = beta1 - 1.96 * se1
qui gen plusse2 = beta2 + 1.96 * se2
qui gen minusse2 = beta2 - 1.96 * se2
gen xaxis2 = xaxis + 0.1
gen xaxis1 = xaxis - 0.1

/* === FIGURE ============================================================*/
	
twoway ///	
	(rcap plusse1 minusse1 xaxis1 ///
		, sort color(blue)) ///			
	(scatter beta1 xaxis1 ///
		, sort lcolor(black) mcolor(blue) lpattern(solid) msymbol(O)) ///	
	(rcap plusse2 minusse2 xaxis2 ///
		, sort color(green)) ///			
	(scatter beta2 xaxis2 ///
		, sort lcolor(black) mcolor(green) lpattern(solid) msymbol(D)) ///		
	, graphregion(color(white))	///
	xtitle("Quintiles of age") ytitle("Treatment effect") ///
	name(figE1b, replace) ///
	xlabel(1 "18-32" 2 "32-39" 3 "39-46" 4 "46-55" 5 "55-90", angle(55)) ///
	legend(on order(2 4) ///
		label(2 "Low income") ///
		label(4 "High income") ///		
		pos(6) col(2) ///
	) ///
	/*legend(off)*/ yline(0) ylabel(-0.04(0.04)0.16, angle(0)) scale(1.3)	
	
graph export "${figures}figE1b.png", as(png) replace
graph export "${figures}figE1b.pdf", as(pdf) replace
